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ABSTRACT 

Context. The literature is rich in analysis and results related to thermally pulsing-asymptotic gi- 
ant branch (TP- AGB) stars, but the problem of the instabilities that arise and cause the divergence 
of models during the late stages of their evolution is rarely addressed. 

Aims. We investigate the physical conditions, causes and consequences of the interruption in the 
^ . calculations of massive AGB stars in the late thermally-pulsing AGB phase. 

5-^ . Methods. We have thoroughly analysed the physical structure of a solar metallicity 8.5Mo star 

^ , and described the physical conditions at the base of the convective envelope (BCE) just prior to 

divergence. 

Results. We find that the local opacity maximum caused by M-shell electrons of Fe-group ele- 
ments lead to the accumulation of an energy excess, to the departure of thermal equilibrium con- 
' ditions at the base of the convective envelope and, eventually, to the divergence of the computed 

models. For the 8.5Mo case we present in this work the divergence occurs when the envelope 
mass is about 2Mq. The remaining envelope masses range between somewhat less than 1 and 



(N 



, more than 2Mq for stars with initial masses between 7 and IOMq and, therefore, our results are 

relevant for the evolution and yields of super-AGB stars. If the envelope is ejected as a conse- 
quence of the instability we are considering, the occurrence of electron-capture supemovae would 
be avoided at solar metallicity. 



X 



5— ( ■ Key words. Stars: AGB and post-AGB ~ evolution - interiors - winds, outflows 

1. Introduction 

One of the most complicated phases of the evolution of low- and intermediate- mass stars is 
the thermally pulsing-asymptotic giant branch (TP-AGB) stage. The physical conditions that are 
sampled by models of this phase are quite extreme and the recurring thermal instability of the 
heUum-burning shell requires codes to be robust on short timescales during the thermal pulse 
cycle. It is common for convergence problems to arise during the study of such stars, although 
rarely are these problems discussed in the recent literature. The failure to converge usually oc- 
curs during the later thermal pulses, when the envelope mass is relatively low due to mass loss 
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through the AGB phase. We are currently using the Monash stellar evolution code MONSTAR 
( Campbell & Lattanziol2008 ) to study the evolution of super asymptotic giant branch (super- AGB) 
stars and have encountered such a problem. 

When the envelope mass is below » 2 M©, the code's solution typically tends to a configuration 
with a negative gas pressure at the bottom of the convective envelope. We find that this convergence 
problem can be delayed by slowly increasing a, the ratio of the convective mixing length to the 
pressure scale height. Our models are originally evolved with a = 1.75. Typically, with a slightly 
higher a, a few more thermal pulses could be modelled until the same convergence problem occurs 
again. This works until the envelope mass is in the region of 0.5 - IMq, when increasing a will 
no longer avoid t he problem. A t ypical value of a ~ 4 is needed to evolve the model up to this 
stage. Previously, iHerwigl (1200 ih has also suggested that a higher a value (of 3) helps to ease 
post- AGB numerical simulations. Recipe s of adopting higher values of a dur i ng the later stag es of 
the AGB evolution are also mentioned in 



Miller Bertolami & AlthausI (l2006h . 



KitsikisI 00081) and 



Weiss & Fergusonl(l2009b . 

In this paper we seek a physical explanation for such behaviour and find that the radiation 
pressure is so large that it supplies all the pressure required for hydrostatic equilibrium. In other 
words, the local l uminosity exceeds the Ed dington limit. This type of instability has previously 
been discussed bv lWood & Faulknen (119861) who suggested that the envelope could be ejected and 
the star would then evolve to the planetary nebula. Eddington limit luminosities as the cause of 
the numerical problems durin g lo w-mass envelope thermal puls es has aheady been identified by 



Lawlor & MacDonaldl (12003 



and 



Miller Bertolami et al 



(l2006l) . 



The stability o f stellar models during the late thermally pulsing AGB phase has als o 



been considered by 


Han et al. 


1994). 


Wagenhuber & Weiss 


(1994 


), 


Wood & Vassiliadis 


( 


1993), 


Renzini et al. 


(1992) 


. Tuchman 


(1984) 


. Barkat & Tuchman 




1980), 


Tuchman et al. 


(1978 


), 


Wood 


J1973L and 


Paczvnski & Ziolkowski ( 


968l). The transition from AGB to post- AGB evolution has 



been considered, from the observational point of view, bv lEngels et alJ (120091) . 

In this paper, we discuss in detail the cause of the instability in the models and whether real 
stars do encounter such an instability in nature. We also discuss possible scenarios for when this 
instability may occur and try to determine whether the envelope would be ejected based on the 
results of our models. The next section of the present work gives a physical description of the 
instability, considering the cases when energy transport is either dominated by radiation or cases 
when convection is efficient. In the third section a hypothesis for the cause of the instability based 
on the existence of an iron element opacity peak is proposed and tested. In the fourth section we 
consider the possible fate of the star after the instabihty. Finally we present the main conclusions 
and discuss our work. 



2. Physical description of tlie instability 

The computation of the late stages of massive TP-(super)AGBs is hampered by the occurrence of 
a type of instability that causes evolutionary models to fail to converge and can affect our under- 
standing of the most advanced phases of the evolution and eventually the resulting white dwarf 



' Note that Table 4 in th is paper is in error, and that all entries for M > 4Mq should be in italics 

3E 



I Lawlor & MacDonak 



20061 : MacDonald, private communication) 
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Fig. 1. Density, temperature, radius, luminosity, k and the ratio of gas pressure to total pressure in 
the instability region. The model profiles represented in solid lines correspond to the times labelled 
in Fig [3] and in particular thick lines represent the last converged models. Dashed lines represent 
the final progress of the instability. 



remnants. This instability is characterised by the development of a zone at the base of the con- 
vective envelope in which the density and the gas pressure, Pgas tend to zero and, therefore, only 
radiation contributes to the total pressure in the equation of state. Eventually the models converge 
to solutions with local Pgas below zero near the base of the convective envelope (BCE) and the 
codes finally crash. As we will show in the next subsection, this is equivalent to having local lumi- 
nosity values higher than the Eddington luminosity. The temperature in the stellar zone in which 
the instability develops is typically 1 - 3 x 10^ K and thus much higher than the values found in the 
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region where hydrogen and helium are partially ionized. Hence the instabihty we are concerned 
with is not directly related to the instability described by IWagenhuber & WeissI (Il994l) . These au- 
thors identified the local decrease in efficiency of the energy transport due to the recombination of 
hydrogen as the main cause of the instability that quenched the computation of AGB evolution. 

We have noticed the presence of the instability described above when using the Mount Stromlo 
and Monash Stellar Evolution code (MONSTA R), but from Siess (private communication) we 
know that the stellar evolution code STAREVOL (ISiessl2010l) . presents the same convergence prob- 
lems during late massive AGB and super- AGB evolution . Both codes and their respective super- 



Dohertv et al 



(120101) . It is worth noting that such convergence 



AGB models have been described by 
problems are common. Many evolution c odes experience numerical pro blems that could be related 



to this instab ility. The MONSTAR code ( Campbell & Lattanzio 



dSiess 



2007 



20101) and the STARS code 



Stancliffe & Jeffervll2007 



I 



2008) and the STAREVOL code 



Lau et al 



20091) have conver- 



gence problems during the late massive AGB or super- AGB evolution. 

We have also checked that whether this instability occurs is insensitive to most input physics 
of the code. For example, we have used different implementations of the OPAL opacities, of the 
treatmen t of mixing (with or without time i nde pendent mixing), or of the mass-loss prescription: 
ReimersI ( ll978h . lVassiliadis & Wood ( 1993 ). or lBloeckerl (Il995br) . The models presented here are 



modelled with the 



Vassiliadis & Wood (119931) mass-loss prescription. These make only minor dif- 



ferences to the occ urrence of the instability. Unless one employs a very high mass-loss rate, such 
as the one used by 



Kovetz et al 



1 2002), the instability cannot be avoided|3 In fact diff'erent input 
physics does change the envelope mass when the instability begins, but they can not prevent it 
from occuring. The opacity calculation could on the other hand, have a huge effect on whether 
this instability occurs, as shown in a later section. Byremoving the Fe opacity peak, it is possible 
to avoid the in stability. Comput ations with old Los Alamos(LAOL) opacities may not suffer this 
instabihty e.g. lBloecken (Il995ah . In fact, the main improvement of the new opacities (OR^.L/OP) 
tables is the treatment of Fe and iron-group opacities in the range of temperatures 10^ - 10^ K 
dlglesias & Rogerslll996al) . This could explain why instabilities are more common with the adop- 
tion of new opacities. This is consistent with our identification of the iron opacity peak as the cause 
of the instability in the models discussed in this paper. 

In terms of numerical details, we have tried varying the obvious things, such as increasing 
the temporal or spatial resolution, but these make little difference. We may delay the instability 
for a few pulses but not eliminate it. However, we discovered that increasing the mixing length 
parameter a enables the code to continue to converge for a lit t le long er. This was also discovered by 



Herwig| ( l200ll) . 



Miller Bertolami & AlthausI (120061) . 



KitsikisI (l2008h and lWeiss & FergusonI (120091) . 



In this way we can push the evolution for a few more thermal pulses, although inevitably the 
convergence problem re-appears. A further increase in a again delays the problem. This is a hint to 
the underlying physics but we are not free to increase a arbitrarily and without limit! 

The instabihty is manifested in the mesh points at the base of the convective envelope. Physical 
conditions at this time are shown in Fig.[T]to|5]for a S.SM© star of solar metallicity. The MONSTAR 
code uses temperature and pressure as dependent variables and the iterations soon produce a T and 
P corresponding to such a high radiation pressure, P^d that the implied gas pressure Pgas - P- Piis 



For their massive stars, 



Kovetz et al 



1 20091) adopted 



Bloecked ( Il995a) mass loss presciption with 77 = 3. 
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Fig. 2. Upper panel: evolution of luminosity associated to H-burning (dotted line) and He-burning 
(solid line); middle panel: evolution of convection; lower panel: evolution of the surface radius, all 
during the second last thermal pulse. Relevant structure profiles of models labelled a to f appear in 
Fig.E] 

negative, as in indicated in Fig.[T] From Fig.|4]we can see that the instabiUty starts at the boundary 
between radiative and convective zones. As the instability develops, the density at the bottom of 
the convective envelope decreases (Fig. [TJ because the model star is expanding, as evidenced by 
the rapid increase in radius a round that region (figure [T]. 



Wood & Faulkner! (119861) . section Illb, while developing an analysis of hydrostatic sequences 
of planetary nebulae, justified the existence of the TP-AGB instability. They computed models 
that would produce the nuclei of the PNe between 0.60 and O.89M0, and realised that, at some 
point in the late evolution of TP-AGB stars, the radiation pressure, Pjad, becomes so large that 
its contribution accounts practically for the total pressure at the base of the convective envelope. 
So Pgas is forced to be zero. These authors suggest that, when this phenomenon takes place in a 
real star, a hydrodynamic expansion of the stellar envelope occurs and that convergence problems 
develop in the models as a consequence of the lack of existence of a hydrostatic solution. 

This is consistent with the instability we are studying in this work. In particular, /3, the ratio 
of gas pressure to total pressure at the base of envelope, drops as the stars evolve along the last 
thermal pulses, with our last successfully converged models yielding values of /3 below 0.001. The 
behaviour of f gas tending to (or falling below) zero is accompanied and, in fact, is equivalent to 
local luminosity values reaching (or exceeding) the Eddington luminosity. We will develop this 
further in the next subsections. 
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Fig. 3. Upper panel: evolution of luminosity associated to H-burning (dotted line) and He-burning 
(solid line); middle panel: evolution of convection; lower panel: evolution of the surface radius, all 
during the last thermal pulse. Relevant structure profiles of models labelled a to f appear in Fig.|9] 



2.1. The importance of the convective efficiency 

The case where the en ergy is carried by radiation alone has already been presented in 



Wood & FauUcnen (119861) The situation is more complicated when convection is present because 
it adds a second energy transport mechanism. The case of vanishingly low convective efficiency 
corresponds to the radiative case. We are interested in how convective efficiency modifies the be- 
haviour of the star. 

In a region where convection is efficient, the local luminosity can exceed the Eddington limit 
without encountering the above instability. This occurs frequently during the lives of stars. The core 
helium flash is one example, where the luminosity can be 10*^ or more, greatly exceeding the 
local Eddington value. Instead convection sets in when the local luminosity exceeds the Eddington 
limit and the highly efficient convection carries most of the energy and the instability does not 
develop. Similarly, during AGB thermal pulses themselves, in the helium burning region we have 
similarly high values and again the efficient convection helps to keep /3 positive even though the 
local luminosity exceeds the Eddington value. 

In a region where convection is efficient, most of the energy is transported by convection instead 
of radiation. We derive below a corresponding Eddington Limit for a convective region and show 
that the value depends on the efficiency of convection. In regions where convection is inefficient 
the local Eddington limit still applies because most energy is transported by radiation. 
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We begin with the equation of pressure equilibrium. 

dP ^ pGM^ 
dr 

where, P is the pressure, r is the radius, G is the gravitational constant and is the mass within 
the radius, and from energy transport by radiation. 
The energy transfer equation is given by 



^ = (2) 

where T is the temperature, k is the opacity, p is the density and is the local luminosity. In the 
radiative case, V = and is defined as • 
The above equation can be rewritten in the form 



_ _ 

dr ~ \6nacr^T^ \Vr 
where 



(3) 



^ _ 3 kL,P ^ 
^ l6nacG mT'^ 



is the gradient if the star is radiative and V is the actual gradient of the star. By definition the 
radiative pressure P^d = ^<^T'^, so we can obtain 

dr W2 ivj 
Combining with ^ = we obtain the diff'erential equation 

dP AncGM, ^ ' 

As we can see from the above equation, a jump in local luminosity or opacity will cause the lo- 
cal radiation pressure to increase much more rapidly than the local pressure. It is therefore possible 
that the radiation pressure completely dominates the total pressure. 

In particular, for the base of the convective envelope (BCE), the stellar region in which we are 
interested, and Mr can be well approximated by the total luminosity of the star, L, and its core 
mass Mc, respectively. If we assume the envelope mass is small and take as the core mass M, 
then the luminosity L is a constant and equals to the value at the base of the convective envelope. 
By integrating over a narrow region around the BCE and using a local average value for the opacity 
< /c>, we can express the ratio of radiation pressure to total pressure at the bottom of the envelope 
as follows: 

£^=i_^.i^m._Li. j_ (7) 

P ^ AncGMc\V,} LEddVr L^^^ 

where <k> = P KdP, which is the average of the opacity over the region. 

From the above equation, we note that yS < is equivalent to L > L!^, where = L^ddY- 
In convective region where convection is efficient V is very nearly equal to the adiabatic gradient 
(Va) and is much smaller than V,, and hence the equivalent Eddington limit is much larger and 
the radiation pressure does not dominate the total pressure. However, in regions where convection 
is inefficient, the value of V approaches the value of Vr, and the above equation approaches the 
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equation ^ - I - /3 ^ 4^cGm ~ I^' Eddington limit converges back to the classical 

value in a radiative region. 

The actual value V depends on the efficiency of the convection and must satisfy the relation 
Vr > V > Va, where Va is the adiabatic gradient. But it is important to keep in mind that radiation 
and convection can coexist in the same region of a star and both can be relatively efficient. It is 
only in some particular regions, such as the stellar centre, in which convection fully dominates and 
rad iation can be ignored. 



Kippenhahn & Weigertl (119941) (page 51) derive the value of V based on mixing length theory. 



The most important quantity is 



3acT^ 8Hp 

where is the mixing length and Hp is the pressure scale-height. The quantity U can be taken as 
the ratio of the radiative 'conductivity' to the convective 'conductivity'. If U is large, convection is 
ineffective and cannot t r anspo rt a substantial fraction of the luminosity. From the equation, also in 



Kippenhahn & WeigertI (11994 



(V-Ve)^/2^ -f/(V,-V) (9) 

where Ve is related to the temperature changes of the moving elements. We see that, in order for 
the equation to have solution when U is larger, V a; and hence most of the energy is transported 
by radiation. 

Because U is inversely proportional to p^, convection is very efficient in the dense central part of 
the star but not necessarily in the envelope where the density is low. In our model, U is of the order 



of 100 at the bottom of the convective envelope, and according to iKippenhahn & WeigertI (11994ft 
this value of U implies radiation dominates over convection for the energy transport. Therefore, 
V ~ V, and the actual Eddington limit is close to the radiative value. When the local luminosity 
exceeds this the instability occurs. 

Notice that this naturally gives us the explanation for why increasing the mixing-length can 
delay the onset of the convergence problems. We see that U is inversely proportional to So by 
increasing a we increase the value of decreasing the value of U, and forcing the convection to be 
more efficient. A smaller U implies more energy is transported by convection. Hence V decreases 
and the effective Eddington luminosity increases beyond the radiative value. Thus by increasing a 
we enable more efficient convection to carry the energy away and the instability is avoided. 

It is interesting to see how the above argument is reflected on an actual computation of a full 
evolutionary sequence. The next few paragraphs compare our previous analysis with the results of 
an actual computation of the final stages of a TP-(super)AGB star. We have considered a S.SM© 
star of solar metallicity, from its main sequence, until the instability is reached after 152 thermal 
pulses. At this time our model star has 2.71Mq, of which I.IIMq correspond to the core and the 
remaining 1.54Mo constitute the envelope. 

Fig. |2] shows the temporal evolution of luminosity, convection and radius during the second 
last thermal pulse. The upper panel shows the temporal evolution of the luminosity associated to 
H-burning, Lu, and He-burning, Lhc, the middle panel shows the evolution of the BCE and the 
convective shell associated with the He-flash and the lower panel shows the variation of the surface 
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Fig. 4. Relevant structure and composition parameters for the model labeled c in Fi^ Upper 
panel: actual gradient (thin solid), radiative gradient (dotted hne) and y8 (thick solid). Middle panel: 
Local luminosity values (dotted line), Eddington luminosity (thin solid) and Eddington luminos- 
ity corrected with the effect of convection (thick solid line). Lower panel: Hydrogen (solid line). 
Helium (dotted line) and Carbon (dashed line) composition profiles. 

radius with time. The labelled vertical lines correspond to selected models whose structure profiles 
are shown in fig |8] and |9] 

Fig.[3]shows the same as the former figure, for the case of the last thermal pulse. As we can see, 
when Liie decays after the flash Lh increases and tends to recover its interpulse value. But, opposite 
to what happens in former thermal-pulses, Lh reaches a maximum of about IQ^^Lq and then de- 
creases again. This is because the stars tries to expand around the region where the gas pressure is 
approaching zero. As a result, the region around the burning shell cools, and the hydrogen burning 
shell cannot reignite. 

Fig. |4] shows a section of the mass profile of the star near the BCE during the last converging 
thermal pulse (at time labelled c in Fig. [3]l. Fig. |5]shows the same as Fig. |4]for the time labelled f, 
that is, for our last converged model. The upper panel of both figures represent the actual gradient, 
V (thin solid line), and the radiative gradient, Vi (dotted line). Convection dominates the energy 
transport in the regions where Vr » V. The thick solid line represents the value of /3 - As we 
can see /3 approaches zero at the bottom of the convective envelope. The middle panels represent 
the local luminosity, (dashed line), the Eddington luminosity computed for radiative regions 
(thin solid line) and the Eddington luminosity computed for convective regions (thick solid line). 



9 



H.H.B.Lau, P.Gil-Pons, C.Doherty, J.Lattanzio: The end of super and massive AGB stars 




1 1-- I I I I L_Lj I I I I J 

1.17 1.1705 1.171 



Fig. 5. Relevant structure and composition parameters for the model labeled / in Fig|3] Upper 
panel: actual gradient (thin solid), radiative gradient (dotted hne) and p (thick solid). Middle panel: 
Local luminosity values (dotted line), Eddington luminosity (thin solid) and Eddington luminos- 
ity corrected with the effect of convection (thick solid line). Lower panel: Hydrogen (solid line). 
Helium (dotted line) and Carbon (dashed line) composition profiles. 

that is, multiplied by the factor ^^Edd.conv -in red. In the convective region, the local luminosity 
could be higher than the Eddington luminosity, but lower than LEdd.conv, as we have justified above. 
However, in the last converged models (see Fig.|5j, ~ i^Edd.conv and y6 ~ in the convective zones 
and we encounter the convergence problem. Note that the instability is actually occuring in a very 
narrow region of about 1 .5 x lO^^M© near the base of the convective envelope. 

Note that both the hydrogen and the heUum burning shells are active at the time of the evolution 
labelled c. But in the very last converged model, labelled f, the hydrogen burning shell is completely 
switched off. Finally, for a reference, the lower panels represent the abundance profile of H, ''He 
and i^c. 

In a summary, our analysis of the convergence limits in our ID -stellar evolutionary code allows 
us to explain a number of consequences that we actually encounter when we perform full evolu- 
tionary computations of super-AGB stars. The code fails to converge when the envelope masses 
are still relatively large (even above 2Mq). The divergence, associated with the increase of lumi- 
nosity in the intershell region above Lndd, occurs in narrow region near the base of the convective 
envelope. Failure of convergence is temporarily solved by increasing the mixing length parameter 
through an increase in a. This increase ultimately allows an increase in the actual Eddington lu- 
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Fig. 6. Left upper panel: mass coordinate versus density near the /fpe peak. Left lower panel: mass 
coordinate versus temperature near the /cpe peak. Right panel. Opacity profile versus mass. The 
insert corresponds to the /fpe peak. 



minosity -conveniently modified by a factor and, therefore, convergence for higher values of 
luminosity in the intershell region. 



3. Test of the Fe-peak opacity hypothesis 



One possible cause for the instability that stops the evolution of our model massive AGB and 
super- AGB stars is the effect of the iron opacity peak (/cpe peak). The A-pe peak is not a unique 
feature in OPAL tables ( Iglesias & RogersI 1996b ). This feature can also be found in opacities from 
the Opacity Project (OP). Seaton & BadnelJ ( 2004 ) has a detailed comparison between OPAL and 
OP opacities. Typically, the Fe-bump occurs at a temperature around 2.0 x 10^ K, this feature being 
more pronounced as the metal content increases. The Fe-bump occurs at a higher temperature in 
OP than in OPAL Jjefferv & Saiolbood) . 



Petrovicetal 



(120061) . was to decrease 



The effect of this /fpe peak in the WR stars, studied by 
the efficiency of energy transport near the local opacity maximum. For stars more massive than 
about I5M0 and close to the Eddington limit such a decrease in the efficiency of energy trans- 
port would lead to local values of the luminosity above Ledd in the zones just below the opacity 
peak. These authors consider that such situation would lead to a significant outwards acceleration 
and the inflation of the stellar envelope. Even though the stars we are considering here are in the 
intermediate-mass range and are in a different evolutionary stage, we also find the /cpe peak (see 
fig|6la nd|7l) at a temperature of fpeak=l -6 x IQr'K, so also close to the value found by 



Petrovic et al 



(l2006h . Furthermore our last converged models also show a fast increase in the surface radius from 



11 



H.H.B.Lau, P.Gil-Pons, C.Doherty, J.Lattanzio: The end of super and massive AGB stars 




Log R/Rg 



Fig. 7. First panel: LogP^ns', second panel: LogT; third panel: k; fourth panel: ebind versus the loga- 
rithm of radius, LogR near the zone of the instability for some selected models prior to the diver- 
gence. The arrows indicate the direction in which time increases. 



1OOO7?0 to about 1700i?o in a time interval of 10 years approximately. Therefore, our model star 



also experiences inflation. iGrafener et al.l (120121) have recently studied inflation in Wolf-Rayet stars 
and luminous blue variables and established a relation with the topology of the /cpe peak. 

In our case the A-pe peak develops after the helium flash and advances inwards, both in mass 
and radius coordinates, as we can see in panel 3 of Fig. |7] In this figure we can see the profiles 
of temperature, gas pressure, opacity and binding energy versus the logarithm of the radius, for 
some selected models between the last helium flash and the time of code divergence. As we have 
explained, the opacity peak is associated to a certain temperature Tpeak- This means that an advance 
inwards of the /cpe peak is directly connected to an advance inwards of Tpeak and, therefore, to the 
cooling and expansion of the afl'ected zone of the star The existence of /cpe peak reduces drasti- 
cally the efficiency of energy transport and thus favours the trapping of the energy generated at the 
helium-burning shell. This implies a local increase of internal energy Mint and, as can be seen in 
panel 4 of Fig. |2]an increase in binding energy ebind, where ebind = fgrav + Mint, and egjav is the gravi- 
tational energy. Mint increases at a zone that is expanding while cooling is caused by recombination, 
and this process is fed-back, that is, more cooling favours more recombination that implies more 
energy for expansion and cooling and so on. 

Eventually the densities between the helium-burning shell and the BCE decrease so much and 
so fast that f gas shows an inversion, that is, an increase with increasing radius that we can see in 
the first panel of Fig. Q This can be identified as the start of the instability that finally leads to the 
crash of the code. The fact that the instability is associated to a decrease in the efficiency in energy 
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Fig. 8. First panel: logarithm of temperature, second panel: specific binding energy ebind, third 
panel: specific internal energy and fourth panel: logarithm of opacity profiles for a few selected 
models along the second last thermal pulse, represented in Fig.|2l 

transport is supported by the fact that the first model presenting f gas inversion corresponds to the 
time at which the inner convective shell disappears. 

This explanation is consistent with section 2. On one hand, the zone of instability identified as 
the zone of reduced density and Pgas -and even Pgas inversion, coincides in mass with the zone of 
instability where [i tends to zero. Furthermore, the peak in binding energy also coincides with the 
zone where > Ledd.conv On the other hand, the instability described in this section can be delayed, 
not only by artificially removing the opacity peak /cpe, but also by increasing the a parameter Both 
increasing a or decreasing k would allow increasing the efficiency of energy transport and the effect 
of blocking of energy would be milder 

The local maximum in the opacity /cpe might be responsible for a behaviour similar to the one 
caused by the /f-mechanism in Cepheids and RR-Lyrae. For these pulsating stars, the hydrogen or 
helium ionisation regions are responsible for a local opacity increase and the loss of efficiency in 
energy transport. This causes an energy excess below the opacity maximum that, eventually, can 
be transformed into work of expansion in the envelope. Once the envelope expands, density and 
temperature decrease and the degree of ionisation also decreases, the star contracts again and the 
process repeats. 

As our code is not hydrodynamic, we cannot follow the e volution furthe r . Neve rtheless it seems 



reasonable to assume that, as in the massive stars studied bv lPetrovic et al.l (12006b , the energy ac- 
cumulated below the /cpe peak will eventually be released and transformed into work of expansion, 
so the (super)-AGB envelopes will also inflate. During this phase of inflation the mass-loss rates 
will be correspondingly high and therefore an enhanced superwind is likely to develop. 
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Fig. 9. First panel: logarithm of temperature, second panel: specific binding energy ebind, third 
panel: specific internal energy and fourth panel: logarithm of opacity profiles for a few selected 
models along the last thermal pulse, represented in Fig. |3] 



We have tested the hypothesis of the /cpe peak being responsible for the divergence of our code 
by artificially removing this peak. We have used fig |6] and Q and inspected the data files of the 
structure profiles of the models prior to divergence. Using this information we have imposed that 
the opacity for temperatures between 1.5 x 10^ K and 3 x 10^ K and densities between 10"" 
and 5 x 10"^g/cm^ must have a constant value equal to the opacity at 3 x 10^ K and density 
5 X 10"^g/cm^. By imposing this constant value, we avoid the instability. This suggests the Fe peak 
is a major reason for this instability, but we must take into account that the A:-peak is real and, unless 
the energy transport is much more eflicient than what we expect from the mixing length theory, the 
peak does exist and is likely to lead, as we have seen, to the inflation and, perhaps, even to the 
disruption of the star 



4. What will happen to the star? 



Wood & Faulknen (Il986h emphasise that the instabihty that affects TP- AGB stars with degenerate 



cores more massive than O.89M0 results from the disappearance of a hydrostatic solution to the 
stellar structure equations. Hence they postulate that the resulting hydrodynamic event will remove 
the stellar envelope. 

We have estimated the outward velocity of the envelope by comparing the radii at the same 
mass mesh points between the last few converging models. We found that almost throughout the 
whole envelope, the shell velocities exceed the local escape velocities. But it is important to note 
that, because the code crashes during the ejection (as seen from the expansion of the radius of the 
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star), the actual velocity during ejection could be even be higher than the one we are estimating. 
This problem, in fact, should be tackled by using hydrodynamical calculations. 

We have also approached the problem of the outcome of the instability by considering the 
binding energy of the envelope. Fig.|8]shows the profiles of Log T, specific binding energy ebind - 
Mint + fgrav, spccific internal energy Mint and opacity LogA: for the models of the second last thermal 
pulse labelled in Fig. |2] and using the same color code. Fig. |9] shows the same for some selected 
models of the last thermal pulse, as labelled in Fig. |3] The specific internal energy Mi„t has been 
computed taking into account the terms asso ciated to the idea l gas, to radiation, to ionization and 



dissociation and to electron degeneracy, as in lHan et al.l (Il994l) . Wherever ebind has a high negative 



value the stellar shells are strongly bound and wherever the binding energy approaches zero, the 
shells are loosely bound. Positive values of the binding energy correspond to the case in which the 
shells would have enough energy to escape, which can be seen in Fig. |7]and Fig. |9] The positive 
values in ebind are located between the top of the helium-burning shell (I.IVOOMq approximately), 
and with the base of the convective envelope (I.I7O6M0 approximately) in Fig.|9] 

The most representative features of the internal energy profiles (third panels of fig|8]and|9]l are 
the the local maxima that correspond to the zones where energy is accumulated due to the /fpe peak, 
as we have explained in the former section. As we can see, these peaks practically coincide in mass 
point with the drop in temperature at the base of the convective envelope and the jump in e^nd- 

The second last thermal pulse is characterised by the fact that the mass point of the base of the 
convective envelope and, therefore, also the Mint peaks and e^nd and k peaks are located practically at 
the same mass point before the helium flash and until the inner convective shell disappears. Along 
the last thermal pulse, though, these masses reach values closer to the centre. Finally, when the 
base of convection reaches the mass point M - \ .\1Q56Mq and temperatures between 0. 1 MK and 
1.5 MK, the opacity peak narrows but reaches a diverging value. Due to this Mim diverges, and so 
does ebind- This is the point at which the it is not possible to proceed further with our calculations. 

Another important feature that distinguishes the second last thermal pulse from the last one 
is the variation of the surface radius (see fig |2]and |3). The second last thermal pulse behaves in 
a standard way. When the helium flash develops the HeBS expand and pushes the HBS that also 
expands and cools down. This leads to the switching off of the HBS and, without its energy supply 
gravity dominates and the stellar envelope contracts. The stellar radius only recovers and increases 
its value when the helium flash is over and the luminosity provided by hydrogen burning is again 
higher than the luminosity due to helium burning. On the other hand, during the last thermal pulse 
the stellar radius remains practically constant. Therefore the energy provided by helium burning 
is not absorbed in the upper HeBS and intershell region. Instead, it is sufficient to sustain the star 
against gravity. This situation is stable as long as the inner convective shell is present. But once this 
disappears while helium burning is still strongly active, a runaway process occurs that leads to the 
increase of the stellar radius on dynamical time scales and, eventually, to the instability described 
in this work and to the divergence of the calculations. 

From this simple approach, we can derive the same conclusion as IWood & Faulknen (Il986l) . 
that is, the envelope will be ejected shortly after the instability is reached. In any case, we would 
like to point out the limitations of our hydrostatic approach to a problem that would be bet- 
ter treated with hydrodynamical codes. Besides, it is also important to realise that the envelope 
masses we are considering in this work -< 2Mq, are much higher than the ones considered in 
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Wood & Faulknen (Il986l) . Therefore, the amount of work necessary to eject completely the stellar 
envelopes of massive AGB stars is much higher and, even though our panels showing the binding 
energy profiles seem to point to the disruption of the star, we cannot discard an eventual fallback 
of the envelopes. We have calculated the envelope binding energy for the last computed models 



as 



Meng et alJ (120081) . and obtained a result of -3x10 erg. This value is negative and, therefore 



technically the envelope still remains bound. But this value should be considered in perspective. 
First, the total binding energy of the star is of the order of 1 0"*^. Therefore the envelope mass (about 
half the mass of the star at this point) has only between 10 and 10"^ of the total binding energy, 
so in practice the envelope is very loosely bound to the core. Second, we can see that the binding 
energy becomes less negative during the last computed models. Thus the envelope becomes less 
bound with time, until the evolution stops and we cannot proceed further, but this trend seems to 
be a good hint of the future behaviour of the star A hydrodynamic code would be more suitable to 
accurately model this ejection process. 

If our estimates showing that most of the envelope (about < 2Mq) will be ejected are correct, the 
outcomes of our model stars might be relatively massive planetary nebulae, with some hydrogen- 
rich matter surrounding a central star of around 1 Mq. However, it is not clear for how long this 
post-ejection system could be observed as planetary n ebula because it could mov e very quickly 
across the HR-diagram to the white dwarf cooling track ( Vassiliadis & Woodlll994 ). The evolution 
of a IMq post- AGB remnant may be too fast to light a massive PN, and even if all massive AGB 
stars undergo the ejection mechanism, they are unlikely to be observed as a PN with a massive 
nebula. Yet, it is interesting to note that the observational counterpart of such a system could indeed 
exist. Planetary nebula N66 contains about 0.6Mq of hydrogen-rich matter i n the nebula, and th e 
lumino sity of the post- AG B star corresponds to a core mass of about 1.2Mq ( Hamann et aklboOS ). 



Indeed, 



VassiliadisI (119961) has already suggested that this planetary nebula could be caused by the 



radiation ejection mechanism. The very low frequency of N66-like PN with a large nebula could 
be consistent with the fact that not all AGB stars going through such an ejection mechanism would 
be observed as a PN 



5. Discussion and conclusions 

In this work we consider the problem of the instability of massive AGB and super- AGB stars and 
identify the reasons why stellar evolution codes break when envelope masses are still relatively 
massive -as massive as 2Mq for the case of super-AGB stars. If this instability is actually en- 
countered by real stars and if it immediately causes their total disruption then it would affect the 
expected yields of intermediate-mass stars. Furthermore, the results of our analysis may have an 
effect on the determination of the mass limit for the formation of electron-capture supernovae, as 
the presence of the instability we have described might quench the evolution of the star toward core 
masses above Chandrasekhar mass. Such an analysis might be particularly important in the case of 
extremely metal-poor super-AGB stars, in which computational time is still a problem and whose 
final fate determination has been often affected by very simple extrapolations. 

The problem of the i nstability of stars at the end of the AGB phase has been considered by 



Wood & Faulknen (119861) . These authors realised that, near the maximum of a late helium flash 



during AGB evolution, local values of the radiation pressure became so high that they were enough 
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to balance gravity. Under these circumstances the gas pressure tended to zero and the local lumi- 
nosity (L,) became higher than the Eddington luminosity (LEdd)- Therefore the model stars departed 
from hydrostatic equihbrium and their evolution was quenched. We have computed the evolution 
of massive AGB stars and realised that the zones at which > Ledd occurred were at the bas e of 
the convective envelope. Therefore we have generalised the results bv lWood & Faulkner! (1 19861) for 



the case in which the departure of hydrostatic equilibrium developed in the presence of convection 
and have generalised an expression of Ledd for convection-dominated environments. 

Opacities ha ve frequently been suggested as re sponsible for the instability. F or instance, in the 
classic works bv lBaker & KippenhahnI (Il965h and lBaker & KippenhahnI (Il965h . it is justified that 
the conditions for the stability of static, radiative layers in gas spheres only depend on the local 
thermodynamical state of these layers -Baker's one-zone model. Thus, if the constitutive relations 
- equations of state and Rosseland mean opacities - are specified, the stability conditions can be 
evaluated without specifying further properties of the layer In the case of solar-composition gas 
spheres the /f-mechanism can even work in regions where H2 dissociation and H-recombination 
are occuring, and such regions are much more extend e d than those of the second ionization zone 
that drives Cepheid pulsations. IWagenhuber & WeissI (119941) also pointed to H-recombination as 
responsible for the termination of the AGB phase. 

We show that, in the case of AGB and super- AGB models at the late stages of their thermally 
pulsing phase a similar instability occurs, but it is caused by the Fe-peak in the opacity at the 



base of the convective envelope. Thi s phenomenon was already described by 



Petrovic et al 



for the case of massive stars. As in 



(I2OO6I) 



Wood & Faulkner! (119861) . the radiation pressure completely 



dominates the gas pressure and density at the convective envelope converges to zero. When this 
instability occurs depends on the treatment of convection. A more efficient scheme of convection 
or overshooting could delay or even avoid its occurrence. In particular, the value of the mixing 
length parameter a determines when the instability occurs. This point is very relevant, as it has re - 
cently been suggested that a should vary during the evolution of the star - iMeakin & ArnettI (l2007h . 
Moreover, i t is possible that mixing length theory cannot adequately describe the convection of the 
actual star - IStanchff^e et alj ( 1201 In . Alternatively, if we could find some observational constraints 
on these models, we might be able to use these to discriminate between different convection theo- 
ries or different values of a for this phase of the evolution. 

We have seen that for massive AGB stars or super-AGB stars with zero-age main-sequence 
masses 7 - IOMq the instabilities occur when the envelope mass is 1 - 2Mq. The instability could 
occur for initial masses as low as 2.5Mq in the solar metallicity case. The stars encountering the 
instability in this study are all hot bottom burning AGB stars. Therefore, the envelope ejected 
should be nitrogen rich. The envelope mass when instability occurs decreases with the masses of 
the AGB. For low mass AGB stars, whether the instabilities occur depend on the choice of the 
mixing length parameter a and treatment of convection. It is fair to conclude the instability does 
occur for massive AGB stars and super-AGB stars unless mixing length theory is fundamentally 
flawed. 

Let us speculate briefly on the behaviour of the stellar material in the region of the instabil- 
ity. Our estimates of the velocity of the convective envelope are greater than the escape velocity. 
Furthermore, when we consider the binding energy profiles for our last converged models we can 
see that binding energy reaches positive values. This suggests that almost all of the envelope could 
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be ejected. In fact the radius of the whole star increases very rapidly. It is thus possible that the 
AGB stars will enter post- AGB phase right after the instability. As to the effects of our analysis 
on the computed yields of AGB and super-AGB stars, we can expect that the ejection phase, if 
it occurs, will be very short compared to the whole TP-(super)AGB phase and that, at this point, 
the nuclear reactions do not play an important role. The nucleosynthesis of the AGB stars will be 
truncated due to the ejection and the yield of the s-process isotopes will be lower compared to the 
alternative assumptions that thermal pulses will keep occuring. 

If ejection does occur after the instability, the TP-(super)AGB evolution would be truncated 
much earlier than in a model star that underwent 'normal' evolution . This would prevent t he oc 
currence o f electron-captured supe movae. Compared to the result of 



Gil-Pons et al. 



(l2007h and 



Poelarends et al 



(l2007h . 



Siess 



(l2008h . which deduce the possibility of electron-captured supernovae 
based on the assumption of regular thermal pulses and extrapolation on core growth rate and mass- 
loss rate, the actual ranges for the occurrence of electron- captured supemovae could be much nar- 
rower in terms of both mass and metallicity. For example, our models show that electron-capture 
supemovae is very unlikely in solar metalhcity. On the other hand, if the instability merely leads 
to a higher mass-loss rate instead of the ejection of the whole envelope, then electron-capture su- 
pemovae may still be possible even though the instability is recurrent. The effect of the instability 
is certainly a factor that cannot be discarded when considering the possibility of electron-captured 
supernovae. On the other hand, at lower metallicities, the opacity jump due to the peak would 
be reduced because of the lower abundance of iron, thus the instability might be reduced and the 
ejection occurs much l ater or not occur at all for lo wer metallicity models. In such a case, if the 
instability described bv lWagenhuber & Weissl(ll994l) cannot operate either, it might possible to find 
a critical maximum metallicity for electron-capture supernovae to occur Evolution of metal-poor 
super-AGB stars will be modeled in future work. 

As a final note, our hydrostatic code could not model through the whole episode of the instabil- 
ity. Although our estimate shows the envelope is Ukely to be ejected, there are a few processes that 
cannot be modelled properly, for instance, the situation when the envelope expands and the density 
converges to zero at the bottom of the convective envelope. As a result, the density of the envelope 
is inversed from the bottom and induce the Rayleigh-Taylor instability, that could facilitate mixing 
of material and energy transport. Thus the instability may be overestimated by our code. However, 
modelling such a process is beyond the cap ability of the Monash evolutionary code. Modelling on 
a 3D hydrodynamic code such as Djeuhty (IDearborn et al.ll2006l) may be needed in order to study 
this instability properly. 
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